import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint
import random

# Generate random numbers between a and b
def GetRandomFactor():
    random.seed()  # Initialize the random number generator with system time, etc.
    return random.uniform(0.9990, 1.0001)

def GetRandomNumber():
    random.seed()  # Initialize the random number generator with system time, etc.
    return random.uniform(-0.1, 0.1)

# Set Chinese font
plt.rcParams['font.sans-serif'] = ['SimHei']  # Use SimHei font
plt.rcParams['axes.unicode_minus'] = False  # Solve the problem of displaying minus sign

# Define seasonal function
def season_factor(t):
    # Divide the year into four seasons, each season is 3 months
    season = (t % 12) // 3  # Each season is 3 months
    return season

# Define time-varying parameter
def r(t):
    season = season_factor(t)

    if season == 1:
        increase_rate = 0.2
    elif season == 2:
        increase_rate = 0.24
    elif season == 3:
        increase_rate = 0.25
    else:
        increase_rate = 0.16
    return increase_rate * 1.1 + 1.2 * np.sin(0.005 * t)

# Define the system of differential equations
def evolution_forests_to_farmland(y0, t, params):
    P, I, B, S, AA, BB = y0  # Set initial population sizes
    K, a, C, a_P, b, I_M, yp_I, rr, yp, Ba, B_M, mu_B, k, S_M, lam, a_AA, m, et_BB, r_AA, K_AA, n, mu_AA, r_BB, K_BB, mu_BB = params

    '''
    K: Plant carrying capacity
    a: Pesticide impact coefficient on plants
    C: Chemical agent concentration
    a_P: Plant resistance coefficient to insects
    b: Insect growth rate
    I_M: Maximum insect population carrying capacity
    yp_I: Insect natural mortality rate
    rr: Chemical impact coefficient on insects
    yp: Bird predation efficiency on insects
    Ba: Bird growth rate
    B_M: Maximum bird population carrying capacity
    mu_B: Bird natural mortality rate
    k: Soil recovery rate
    S_M: Maximum soil fertility
    lam: Chemical impact on soil
    a_AA: Species AA grazing coefficient on plants
    m: Species AA competition coefficient with insects
    et_BB: Species BB predation coefficient on insects
    r_AA: Species AA growth rate
    K_AA: Species AA carrying capacity
    n: Insect competition coefficient with species AA
    mu_AA: Species AA natural mortality rate
    r_BB: Species BB growth rate
    K_BB: Species BB carrying capacity
    mu_BB: Species BB natural mortality rate
    '''
    dPdt = r(t) * P * (1 - P / K) - a * P * C - a_P * P * I - a_AA * P  # Plant growth
    dIdt = b * I * (1 - I / I_M - m * AA / I_M) - rr * I * C - yp * B * I - yp_I * I - et_BB * I + (P / 300) * 15  # Insect growth
    dBdt = Ba * B * (1 - B / B_M) + yp * B * I - mu_B * B  # Bird growth
    dSdt = k * (S_M - S) - lam * S * C  # Soil quality
    dAAdt = r_AA * AA * (1 - AA / K_AA - n * I / K_AA) + a_AA * P - mu_AA * AA  # Species AA growth
    dBBdt = r_BB * BB * (1 - BB / K_BB) + et_BB * I - mu_BB * BB  # Species BB growth
    return [dPdt, dIdt, dBdt, dSdt, dAAdt, dBBdt]

# Initial population sizes
P0 = 350  # Plants
I0 = 130  # Insects
B0 = 30  # Birds
S0 = 1.2  # Soil fertility
AA0 = 37  # Species AA
BB0 = 42  # Species BB
y0 = [P0, I0, B0, S0, AA0, BB0]

# Time points
t = np.linspace(0, 150, 1000)

# Parameter design
params = (
    500,  # K: Plant carrying capacity
    2.8,  # a: Pesticide impact coefficient on plants
    0.01,  # C: Chemical agent concentration
    0.0005,  # a_P: Plant resistance coefficient to insects
    0.45,  # b: Insect growth rate
    200,  # I_M: Maximum insect population carrying capacity
    0.05,  # yp_I: Insect natural mortality rate
    20.8,  # rr: Chemical impact coefficient on insects
    0.0007,  # yp: Bird predation efficiency on insects
    0.13,  # Ba: Bird growth rate
    35,  # B_M: Maximum bird population carrying capacity
    0.08,  # mu_B: Bird natural mortality rate
    0.13,  # k: Soil recovery rate
    1.5,  # S_M: Maximum soil fertility
    2.8,  # lam: Chemical impact on soil

    0.015,  # a_AA: Species AA grazing coefficient on plants
    0.01,  # m: Species AA competition coefficient with insects
    0.15,  # et_BB: Species BB predation coefficient on insects
    0.2,  # r_AA: Species AA growth rate
    50,  # K_AA: Species AA carrying capacity
    0.02,  # n: Insect competition coefficient with species AA
    0.18,  # mu_AA: Species AA natural mortality rate
    0.3,  # r_BB: Species BB growth rate
    50,  # K_BB: Species BB carrying capacity
    0.3  # mu_BB: Species BB natural mortality rate
)

# Solve the system of differential equations
solution = odeint(evolution_forests_to_farmland, y0, t, args=(params,))

# Plot the results
plt.figure(figsize=(10, 8))
plt.plot(t, solution[:, 0] / P0, label='Current biomass of plants / Initial biomass (P)')
plt.plot(t, (solution[:, 1] / I0) * GetRandomFactor() + GetRandomNumber() * 0.55 * GetRandomFactor() * 1.2 * np.sin(t), label='Current biomass of insects / Initial biomass (I)')
plt.plot(t, solution[:, 2] / B0 + GetRandomNumber() * 0.05 * GetRandomFactor() * np.sin(t), label='Current biomass of birds / Initial biomass (B)')
plt.plot(t, solution[:, 3] / S0, label='Soil fertility / Initial soil fertility (S)')
plt.plot(t, solution[:, 4] / AA0 + GetRandomNumber() * 0.05 * GetRandomFactor() * np.cos(t), label='Current biomass of herbivorous species A / Initial biomass (A)')
plt.plot(t, solution[:, 5] / BB0 + GetRandomNumber() * 0.05 * GetRandomFactor() * np.sin(t), label='Current biomass of bird species B / Initial biomass (B)')

plt.legend()

plt.xlabel('Time (months)')
plt.ylabel('Initial biomass (soil fertility) / Current biomass (soil fertility)')
plt.title('Population size/soil fertility over time after farmland ecosystem maturation (using chemicals)')
plt.show()